Collapse and revival of oscillations in a parametrically excited Bose-Einstein 
condensate in combined harmonic and optical lattice trap 
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In this work, we study parametric resonances in an elongated cigar-shaped BEC in a combined 
harmonic trap and a time dependent optical lattice by using numerical and analytical techniques. 
We show that there exists a relative competition between the harmonic trap which tries to spatially 
localize the BEC and the time varying optical lattice which tries to delocalize the BEC. This 
competition gives rise to parametric resonances (collapse and revival of the oscillations of the BEC 
width). Parametric resonances disappear when one of the competing factors i.e strength of harmonic 
trap or the strength of optical lattice dominates. Parametric instabilities (exponential growth of 
Bogoliubov modes) arise for large variations in the strength of the optical lattice. 
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<*2 . I. INTRODUCTION 

A periodic variation of a system parameter leads to the exponential growth of certain modes of the system. This 
phenomenon is known as parametric resonance and occurs in a wide variety of systems such as classical oscillators, 
nonlinear optics, system governed by non-linear Schrodinger equations and in Hamiltonian chaotic systems. Para- 
metric excitations in the form of Faraday waves, were observed in a cigar-shaped Bose-Einstein condensate (BEC) by 
periodically modulating the radial trap frequency [|| ■ The modulation of the radial trap frequency leads to a periodic 
modulation of the density of the cloud in time, which in turn leads to a periodic change in the nonlinear interactions. 
This leads to the parametric excitation of longitudinal sound-like waves (Faraday waves) in the direction of weak 
confinement. From the theoretical side, parametric excitations have already been investigated by numerous authors 
*^ \ [341) HHUI ■ A BEC in an optical lattice potential exhibits parametric excitations when the depth of the optical lattice 
Ch ■ potential is periodically modulated in time. These parametric excitations correspond to an exponential growth of the 
population of counter-propagating Bogoliubov excitations [HI, [l3[ . When the frequency of oscillation of the scattering 
i ^ i , length is an even multiple of the radial or axial natural oscillation frequency, respectively, the radial or axial oscillation 
of the condensate exhibits resonance 5]. Parametric resonances in a BEC under periodic oscillations of the position 
of a quadratic plus quartic trap have also been predicted Q. Experimentally, Stoferle et al. 0, [H| measured the 
^ 1 excitation spectrum of the condensate in the superfluid regime by time modulating the optical lattice depth. The 
width of the expanding cloud after it is released from the trap is taken as a measure of the excitation energy. The 
observed resonances in these experiments were explained by the generation of parametric amplification of Bo goli ubov 
states and the subsequent nonlinear dynamics, which leads to the broadening of the momentum distribution |12|. [l3|. 
• , In this work, we study parametric resonances in an elongated cigar-shaped BEC in a combined harmonic trap and 
a time dependent optical lattice by using numerical and analytical techniques. In order to distinguish the present 
work from parametric resonances studied earlier, we always consider the case when the frequency of oscillation of the 
optical lattice depth is not equal to the axial or radial natural oscillation frequency. We show that there exists a 
relative competition between the harmonic trap which tries to spatially localize the BEC and the time varying optical 
lattice which tries to delocalize the BEC. This competition gives rise to parametric resonances (collapse and revival in 
the oscillations of the BEC width). Parametric resonances disappear when one of the competing factors i.e strength of 
■ harmonic trap or the strength of optical lattice dominates. On one hand, parametric instabilities (exponential growth 
of Bogoliubov modes) arise for large variations in the strength of the optical lattice and on the other hand, stable 
oscillations of the condensate width occurs when the harmonic trap strength dominates. We employ a variational 
technique [l6l - fl8j to derive analytical predictions for the parametric resonances. The analytical results qualitatively 
agree with the numerical results. 



II. MODEL 



We consider an elongated cigar shaped Bose-Einstein condensate (BEC) of N two- level 87 Kb atoms in the \F = 1 > 
state with mass m and frequency uj a of the \F = 1 >— s- \F' = 2 > transition of the line of 87 Rb in combined 
harmonic and optical lattice potential. In order to create an elongated BEC, the frequency of the harmonic trap 
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along the transverse direction should be much larger than one in the axial (along the direction of the optical lattice) 
direction. This system is described by the standard mean- field model for the BEC, the GP equation in 3D which is 
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where tp(r,t) is the condensate wavefunction, m is the mass of an atom, N represents the total number of the 
atoms in the condensate, a is the atomic scattering length which can be tuned by the Feshbach resonance technique, 
g = Ai:h 2 a/m is the three dimensional inter-atomic interaction. The normalization condition of the order parameter 
is J \tp(r,t)\ 2 dr — 1. The confining potential for the BEC is given by 
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where uj 7 
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ku, u) z — Xlo represents the trapping frequencies in the respective directions. We have assumed 



that harmonic potential in the x-direction is superimposed by a ID OL potential of the form 



,2nx 



= Vsin 2 ( — )(1 + asin(nt)) 



(3) 



where V is the amplitude of the OL potential. K = ^ and A; is the wavelength of the laser beam used to create 
ID OL potential. 7 is the parameter representing the strength of harmonic confinement such that < 7 < 1. In our 
all calculations 7 = 0.001 shows very weak and 7 = 1 shows very strong harmonic confinement. f2 is the frequency 
of external perturbation at any time t in the time dependent ID OL potential with the parameter a < 1. We are 
interested in studying the effect of time dependent potential on the dynamics of the cigar-shaped BEC. Taking the 
y,z directions to be the transverse and x as longitudinal we assume the confinement in the transverse directions to 
be sufficiently strong (u z ,u y » u) x ) so that the transverse degrees of freedom (y, z) freeze out and the condensate 
motion gets confined to the longitudinal direction only and the BEC is quasi ID. 

To solve the above problem numerically, we first make the EqnQJlimensionless through a set of linear transformations 

t — > ~ , x — » XI , y — >Yl, z — > Zl , I — J jj^jj , tp(x, y, z, t) — ^ X ^^' T ^ which we will use throughout the paper. 
After these transformations Eqn. Q] becomes 
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where 773D is the non-linearity in the GP equation. To reduce the above Eqn to quasi-lD form it is assumed that the 
system remains confined to the ground state in the transverse direction. By integrating the above Eqn with respect 
to the ground state transverse wavefunctions using the method of separation of variables, setting i/=l,we obtain quasi 
ID form as , 
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where rjio = 2aN^f\~K,/l with the normalization f,°° <fi(X, r)dX = 1 and with v=\ we assume I to be the harmonic 



oscillator's length in longitudinal direction. Also 
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Therefore we can rewrite the above EqnJS] as 
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Figure 1: Parametric resonances and instability in the width of the BEC vs dimensionless time r in a combined harmonic trap 
and a time modulated optical lattice. The parameters are, 7 = 0.02, (al): a = 0.01, (a2): a = 0.1, (a3): a = 0.2,(a4): a — 0.5. 

III. NUMERICAL RESULTS 

In the numerical simulations, we have used the split-step Crank Nicolsan method to solve the GP equation in one 
space variable. We have used the imagetimclD.F and realtimelD.F codes [19j. The ground state solution is found 
by imagetimelD.F code with NPAS = 100000, NRUN = 20000 with space step dX = 0.01, time step dt = 0.0002. 
To find the ground state solution, we have used the initial Gaussian wavefunction 4>{X, r = 0) = exp~ x / 2 /^/n in the 
imagetinelD.F code. Then to study the time evolution of the GP equation, we have used the realtimelD.F code with 
dr = 0.001. The results of the numerical simulations are now presented below. 

To perform a systematic numerical study of the above Eqn. [51 we take the value of A/ = 4.0, V Q = 10.0, uj r = 0.5 
. To perform the simulations we have first prepared the ground state wavefunction of BEC using the imaginary time 
propagation method with the V(X, r) at r = 0. At r = 0, T = 1 and then this ground state is propagated in the 

real time with V{X,t) = ^ + AT sin 2 [^^j. We have used the value of dX = 0.01, dr = 0.0001, rj 1D = 17.0. 

When the frequency of oscillation of the scattering length is an even multiple of the radial or axial natural oscillation 
frequency, respectively, the radial or axial oscillation of the condensate exhibits resonance § . In order to distinguish 
the present work from parametric resonances studied earlier [Bj], we always consider the case when the frequency of 
oscillation of the optical lattice depth is not equal to the axial or radial natural oscillation frequency. The natural 
frequency of axial oscillations is 2lu, so in order to stay away from resonance we have taken uj r = 0.5 as mentioned 
above. 

Fig.l shows the time evolution of the rms value of the condensate axial width < X rms > in a weak harmonic 
potential 7 = 0.02 for four different values of the perturbation strength a — 0.01 (plot al), a = 0.1 (plot a2), 
a = 0.2(plot a3), a = 0.5 (plot a4). We clearly see the competing effect of the harmonic potential and the periodic 
modulation of the intensity of the optical lattice. The harmonic potential tries to localize the condensate at the 
center of the trap while the periodic modulation of the lattice generates Bogoliubov excitations which with time 
grows and tends to delocalize the condensate. In plot al, the periodically varying part of the lattice asinw r r is 
weak and varies between zero and a. The collapse and revival seen in plot al is a result of competition between the 
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Figure 2: Parametric resonances and instability in the width of the BEC vs dimensionless time r in a combined harmonic trap 
and a time modulated optical lattice. The parameters are, 7 = 0.1, (bl): a — 0.01, (b2): a — 0.1, (b3): a = 0.2,(b4): a = 0.5. 



harmonic potential and a sincj r r. Due to the weak value of the perturbing factor a, the harmonic potential does not 
allow the condensate to become unstable and the balance between both the competing factors induces a continuous 
growth and decay cycles in the amplitude of rms value of the condensate width. Now as we increase the perturbing 
strength a to 0.1(a2), parametric amplification of the Bogoliubov modes begin to appear after a certain time interval. 
After a certain time, the harmonic potential is no longer able to stabilize and localize the condensate and there is an 
unprecedented growth of the parametric excitations and the condensate becomes unstable. On increasing a further 
to 0.2 (a3) and 0.5 (a4), the parametric amplification and hence instability of the condensate is more pronounced 
and the growth of the unstable modes starts almost immediately and the harmonic potential is completely unable to 
stabilize the condensate. 

In Fig. 2, we now show the dynamics of the rms value of the condensate width in slightly stronger harmonic 
potential 7 = 0.1 for same values of a, a = 0.01 (plot bl), a = 0.1 (plot b2), a — 0.2(plot b3), a = 0.5 (plot b4). 
Now since in Fig. 2, the harmonic potential is stronger than that in Fig.l, naturally the value of a needed to make 
the condensate unstable and amplify the unstable modes is higher. At a = 0.1, we observe that the width still shows 
the typical growth and decay cycles as compared to the previous case of Fig.l, where at a = 0.1, the condensate 
becomes parametrically unstable after a certain time. In this case now, the condensate is completely unstable only 
when a > 0.2. Another important point to note is that the rms width of the condensate in the relatively stronger 
harmonic potential of Fig. 2 is much smaller as compared to that in Fig.l. 

The parametric instability observed in Fig.l and Fig. 2 for large value of a is a consequence of the exponential 
growth of the Bogoliubov modes of the system induced by the periodic variation of the optical lattice depth. In order 
for these Bogoliubov modes to be amplified, they must preexist at t = 0. This means that the initial state is not a 
pure ground state. Indeed the system is certainly not in its ground state at t — due to excitations produced when 
loading the condensate into the lattice. Some numerical noise is always present and yields a very small occupation of 
Bogoliubov modes. The exponential growth of the parametric instability with time is due to the nonlinear mean-field 
dynamics of the condensate governed by the GP equation. This nonlinear dynamics is triggered by the growth of the 
parametrically unstable modes. 

The growth of the rms width with time due to the periodic variation of the optical lattice also motivated us to 
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Figure 3: Density distribution |0(X)| 2 at r = 1000 vs X corresponding to the cases considered in Fig.l. The parameters are 
7 = 0.001, L(al): a = 0.01, L(a2): a = 0.1, L(a3): a = 0.2, L(o4): a = 0.5. 



analyze how this periodic modulation influences the spatial extension of the condensate density. Fig. 3 and Fig. 4 show 
the density plots corresponding to Fig.l and Fig. 2 respectively. Clearly as the strength of the perturbation a increases, 
the wavefunction also gets progressively delocalized and starts extending over more lattice sites. On increasing the 
strength of the harmonic potential (Fig. 4), the influence of a in delocalizing the condensate wavefunction is less 
effective, since a stronger harmonic potential now localizes the condensate to a greater extent and acts against the 
perturbation. The periodic modulation redistributes the atoms in the lattice wells in a random manner and hence 
the nonuniform density distribution across the lattice wells for large a is seen in Fig. 3 and Fig. 4. 



IV. VARIATIONAL ANALYSIS 



All the results that we have presented in the previous section are fully numerical. In this section, using the Lagrange 
variational technique, we reproduce the essential physics presented in the previous section. 

The dynamics of a trapped one-dimensional repulsive Bose gas is described in the framework of the ID Gross- 
Pitaeviskii equation 
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with the total number of atoms N = J \ip\ 2 dx. Eq. [7]is obtained for a highly elongated cigar shaped BEC. In order 
to create a cigar shaped elongated BEC, the frequency of the harmonic trap along the transverse direction should 
be made larger than the one in the axial (along the direction of the optical lattice) direction. The condition for ID 
approximation is u) y ,u) y uj X) /i <C hui y , hu> z , where \x is the chemical potential. To analytically solve Eq. we need 
a trial wavefunction. A proper choice of the trial wavefunction is very important in this technique, so we are assuming 
that the total wavefunction is a gaussian envelope multiplied by a periodic part. The periodic part of the wavefunction 
which arises from the optical lattice consists of a unmodulated and a modulated part. The modulated part of the 
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Figure 4: Density distribution |0(X)| 2 at r = 1000 vs X corresponding to the cases considered in Fig.l. The parameters are 
7 = 0.1, L(fel): a = 0.01, L(b2): a = 0.1, L(b3): a = 0.2, L{b4): a = 0.5. 



periodic wavefunction is small compared to the unmodulated part because the optical lattice modulates the gaussian 
wavepacket along the longitudinal direction only slightly. To this end, we use the following trial wavefunction for the 
periodic part of the wavefunction 



1 + (72 T' COS 



/47TX\ 



(8) 



where A' — V and T' — 1 + a sin(f2t) and 02 is a variational parameter. 
Using the above ansatz Eq. [8j we can write the trial wavefunction as 
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where W x is the width of the condensate, a x is the slope, f} x is the (square root of the curvature radius) 1 , xq is 
the center of the cloud and <f> is the global phase are the variational parameters. Above trial wavefunction can be 

normalized with respect to N in the limit e ^ — > 0, where j is any positive integer. We are taking this approximation 
in all the calculations because A is very small as compared to W x . This leads us to the final trial wavefunction 
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The corresponding ID Lagrangian density is 
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Figure 5: Analytical results of the BEC width (plots bbl and bb2) vs dimensionless time r and density \ip(X, r)| 2 (L(bbl) and 
L(bb2) at r = 300 vs X. The parameters are a = 0.2, 7 = 0.1 (plots bbl and L(bbl)) and 7 = 0.008 (plots bb2 and L(bb2)). 
Here ^(X,r) = ^^. 



The ID Lagrangian of the system is found from as 

Lid =< £ >= / £ dx. 
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This yields, 
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where A' ,T' are the same as already described and 
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The evolution of the Lagrange parameters, qi = {W x , x a , a x , j3 Xl a 2) is ruled by the corresponding set of Lagrange 
equations 
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which give us equation of motion for the conservation of the norm, 



dt 

the equation of motion for the slope gives, 
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the equation of motion for the center of cloud gives, 



and for width gives, 
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The equation of motion for the parameter <ii yields, 

A'/2 



<J2 = 



V^Wx \ A( J Am 



(22) 



In calculating the above equation of motion we have assumed that <r 2 << 1 so that we have ignored all higher order 
terms. 

On rearranging the above equations we get evolution of the slope and curvature as 

to W T , . 
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and finally the evolution of width as 
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where t]id = ^pr- Let us introduce the constant P = \J~^Nj where P denotes the strength of interaction between 



two atoms. After making the transformations t — > — , x ^ XI , y — > Yl , z —± Zl , I = \ — , the dimensionless 
equation of motion for the width (fx leads to 



W X+1 W X = ^ + 
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where a± = J in the limit uj± » uj x , uj± = uj v ,uj z . Note that lj = u x since v = 1. 
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While solving the differential equation of width we have taken P = 0.07, N = 460, rjm = 17.0 where a — 5.lnm for 
the rubidium atoms , loj_ = IOOcj^,. In Fig. 6, we show the time dynamics of the width calculated from Eqn. [25] in a 
harmonic potential of strength a — 0.2 for two values of 7 = 0.1 (left plot) and 7 = 0.008 (right plot). Qualitatively 
we see the similarity with our numerical results. For strong value of 7 = 0.1, we clearly see the collapse and revival 
in the oscillations of the width while for a weaker value of 7 = 0.008, the parametric instabilities start to grow and 
the width of the condensate increases very fast with time. The density plots corresponding to bbl and bb2 are shown 
in L(bbl) and L(bb2) respectively. Clearly we see that the wavefunction for the case 7 = 0.008 is more delocalized 
as compared to the case for 7 = 0.1. This result is again in agreement with numerical results. The occurrence of 
this type of parametric instability associated with Bogoliubov excitation in Bose- Einstein condensates originates from 
their supernuid nature and the existence of undamped collective excitations. 

V. CONCLUSION 

In this work we have explored the parametric excitations in an elongated Bose-Einstein condensate subject to a 
combined harmonic potential and ID optical lattice whose intensity is periodically modulated in time. On one hand 
the periodic modulation generates Bogoliubov excitations which grows with time and tends to increase the condensate 
size and hence delocalize it and on the other hand the harmonic trap tries to localize the condensate at the center of 
the trap and act against the instability of the Bogoliubov modes. The ensuing nonlinear dynamics leads to parametric 
resonances (collapse and revival of the oscillations of the condensate width) . The parametric resonances and instability 
are nicely accounted for by the numerical and analytical calculations. 
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